% [light, normal] = getLightDirs(grayChrome, grayChromeMask)
%
% input:
%   * grayChrome - The gray image of the chrome with highlight point
%   * represented in a numRows-by-numColumns matrix.
%   * grayChromeMask - The gray image of the chrome mask represented in a
%   numRows-by-numColumns matrix indicating whether a corresponding pixel
%   in the grayChrome image is a scene point or not. However, this matrix
%   doesn't consist of only 0 and 1 indicating exactly whether a pixel is
%   or is not a scene point. The values kind of range from 0 to 255. In
%   this method, if grayChromeMask(i,j) is larger than 127.5, then
%   pixel(i,j) in the image is considered to be a scene point; otherwise,
%   the background pixel.
%
% output:
%   * light - A unit vector indicating the light direction of the image.
%   The vector is of the form (l_x, l_y, l_z), where the x axis is along
%   the column index increasing direction (towards right) and the y axis is
%   along the row index increasing direction (downwards). The z axis points
%   out of the image.
%   * normal - A unit vector indicating the normal direction of the
%   highlight point. The vector is in the form of (n_x, n_y, n_z), where
%   the x axis is along the column index increasing direction (towards
%   right) and the y axis is along the row index increasing direction
%   (downwards). The z axis points out of the image.
%
% description:
% This method takes a gray chrome image with highlight point and a chrome
% image mask to decide the lighting direction of the image.

function [light, normal] = getLightDirs(grayChrome, grayChromeMask)

% find the column index (x) and row index (y) of the highlight point
[yh, xh] = findHighlightPosition(grayChrome);

% find the radius of the chrome ball
[yc, xc, R] = findChromeFeatures(grayChromeMask);

% calculate the lighting direction
[light, normal] = findDir(xh, yh, xc, yc, R);
end

function [r, c] = findHighlightPosition(gray_img)
maxI = max(max(gray_img));
maxPxls = (gray_img == maxI);
[numRs, numCs] = size(gray_img);
rIndices = 1:numRs;
cIndices = 1:numCs;
maxPxls_ri = maxPxls.*repmat(rIndices',1,numCs); 
maxPxls_ci = maxPxls.*repmat(cIndices,numRs,1);
numMaxs = sum(sum(maxPxls));
r = round(sum(sum(maxPxls_ri))/numMaxs);
c = round(sum(sum(maxPxls_ci))/numMaxs);
end

function [rc,cc,R] = findChromeFeatures(gray_mask)
threshold = max(max(gray_mask))/2;
realPxls = (gray_mask >= threshold);
[numRs, numCs] = size(gray_mask);
rIndices = 1:numRs;
cIndices = 1:numCs;
realPxls_ri = realPxls.*repmat(rIndices',1,numCs);
realPxls_ci = realPxls.*repmat(cIndices,numRs,1);
numReals = sum(sum(realPxls));
rc = round(sum(sum(realPxls_ri))/numReals);
cc = round(sum(sum(realPxls_ci))/numReals);
R1 = max(max(realPxls_ri))-rc;
R2 = max(max(realPxls_ci))-cc;
R = round((R1+R2)/2);
end

function [l, n] = findDir(xh, yh, xc, yc, R)
zh = sqrt(R^2 - (xh-xc)^2 - (yh-yc)^2);
n = [xh-xc, yh-yc, zh];   % the normal vector of the highlight point
zr = n*n'/zh/2;
r = [0, 0, zr];     % the reflection vector
l = r - n;
n = n/(sqrt(n*n'));
l = -l/(sqrt(l*l'));
end

